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ABSTRACT 


The  problem  of  loss  of  stability  of  marine  vehicles  under  cross  track  error 
control  in  the  presence  of  mathematical  versus  actual  system  mismatch  is  analyzed. 
For  demonstration  purposes,  variations  in  the  heading  angle  control  gain  are  studied. 
Particular  emphasis  is  placed  on  analyzing  the  effects  of  observer  design  on  system 
response  after  initial  loss  of  stabile  of  straight  line  motion.  It  is  shown  that  the 
dynamics  of  the  observer  may  h-  a.  «v  -  ificant  effect  on  the  computed  gain  margin 
of  the  control  system  depending  on  particular  basis  used. 
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I.  INTRODUCTION 


Accurate  path  control  of  surface  ships  and  underwater  vehicles  along  prescribed 
geographical  paths  is  a  basic  problem  that  is  becoming  increasingly  important,  par¬ 
ticularly  as  the  missions  of  ocean  vehicles  become  more  complicated  with  strict 
requirements  for  performance.  In  order  for  a  control  law  to  be  able  to  perform  its 
mission  in  a  realistic  operational  scenario  it  has  to  be  robust  enough  so  that  it  can 
maintain  stability  and  accuracy  of  operations  in  the  presence  of  modeling  errors  and 
environmental  uncertainties.  The  robustness  properties  of  the  design  are  particu¬ 
larly  important  due  to  the  unpredictable  nature  of  the  ocean  environment  and  the 
changes  in  the  hydrodynamic  characteristics  of  the  vehicle  during  turning,  changes 
in  the  forward  speed,  or  operations  in  proximity  to  other  objects  in  the  area.  For 
these  reasons,  there  exists  a  need  for  the  analysis  of  the  robustness  characteristics 
of  a  particular  control  law  design  and  the  establishment  of  a  rational  operational 
envelope  based  on  stability  and  performance  criteria.  Previous  studies  [Parsons 
and  Cuong  (1977)]  showed  that  gain  adaption  is  highly  desirable  due  to  changes 
in  the  linearized  vehicle  hydrodynamics  with  different  operation  conditions,  such 
as  depth  under  keel.  The  resulting  adaptation  scheme  [Parsons  and  Cuong  (1980)] 
required  significant  vehicle  motion,  which  may  be  undesirable  when  operating  in 
restricted  waters,  or  in  object  recognition  and  localization  tasks.  Integral  control 
techniques  [Parsons  and  Cuong  (1981)]  proved  quite  effective,  but  neglected  the 
nonlinear  behavior  of  the  vehicle,  which  becomes  very  important  at  low  speeds  and 
hover  operations.  Model  based  compensators  exhibit  robust  behavior  under  con¬ 
ditions  of  parameter  uncertainty,  which  is  as  good  as  the  classical  linear  quadratic 
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regulators  for  linear  output  feedback  systems  [Healey  (1992)].  Alternatively,  sliding 
mode  controllers  exhibit  very  robust  characteristics  given  an  estimate  of  the  param¬ 
eter  uncertainty  and/or  disturbances  [Papoulias  and  Healey  (1992)],  [Yoerger  and 
Slotine  (1985)].  Sliding  mode  control,  however,  does  not  offer  an  infinitely  robust 
design  and  it  suffers  from  a  series  of  bifurcation  phenomena  and  loss  of  stability 
unless  proper  care  is  exercised  [Papoulias  (1991)]. 

In  this  work  we  analyze  the  problem  of  loss  of  stability  of  a  marine  vehicle 
under  cross  track  error  control  in  the  presence  of  mathematical  versus  actual  system 
mismatch.  For  demonstration  purposes,  variations  in  the  heading  angle  control 
gain  are  studied.  Previous  studies  [Oral  (1993)]  concentrated  on  system  response 
assuming  perfect  and  complete  state  measurement.  Particular  emphasis  in  this  work 
is  placed  on  analyzing  the  effects  of  observer  design  on  system  response  after  initial 
loss  of  stability  of  straight  line  motion.  The  main  loss  of  stability  cases  analyzed 
here  occur  in  the  form  of  generic  bifurcations  to  periodic  solutions  [Guckenheimer 
and  Holmes  (1983)].  We  use  center  manifold  reduction  techniques  and  averaging 
in  order  to  capture  the  stability  properties  of  the  resulting  limit  cycles  [Chow  and 
Mallet-Paret  (1977)].  It  is  shown  that  the  dynamics  of  the  observer  may  have  a 
significant  effect  on  the  computed  gain  margin  of  the  control  system  depending 
on  the  particular  basis  used.  All  computations  in  this  work  are  conducted  for  the 
NPS  autonomous  underwater  vehicle  [Bahrke  (1992)]  and  all  results  are  presented  in 
standard  dimensionless  quantities  with  respect  to  vehicle  length,  7.3  ft,  and  nominal 
forward  speed,  2  ft/sec. 
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II.  PROBLEM  FORMULATION 


A.  EQUATIONS  OF  MOTION 

The  linear  maneuvering  equations  of  motion  of  a  marine  vehicle  in  the  hori¬ 
zontal  plane  are  written  in  dimensionless  form  as, 

m( v  +  r  +  xGr)  =  Y^r  +  Y„v  +  Yrr  +  Y„v  -|-  YgS  (2.1) 

Itr  -f  mxa( v  +  r)  =  N+r  +  N^v  +  Nrr  +  N„i /  +  NgS  (2.2) 

where  all  symbols  are  explained  in  the  nomenclature.  Equations  (2.1)  and  (2.2)  can 
be  used  to  derive  a  second  order  transfer  function  between  the  rudder  angle  6  and 
yaw  rate  r.  For  low  frequency  maneuvering  motions  this  second  order  equation  can 
be  approximated  by  expanding  in  Taylor  series  and  keeping  the  first  order  terms 
only.  The  result  is 

r  =  or  -f  bS  (2.3) 

Equation  (2.3),  which  is  sometimes  referred  to  as  Nomoto’s  first  order  model,  is 
particularly  useful  in  control  system  design  since  no  sway  velocity  feedback  is  neces¬ 
sary.  This  equation  predicts  linear  variation  of  the  steady  state  turning  rate  versus 
rudder  angle.  In  reality,  the  r-6  curve  has  characteristics  of  softening  spring  mainly 
due  to  speed  loss  during  turning.  To  account  for  this  a  modified  version  of  equation 
(2.3)  is  used, 

r  =  ar  +  a3r3  +  b6  (2.4) 

where  a3  is  usually  determined  from  steady  state  results.  Finally,  the  model  is 
complete  by  the  incorporation  of  the  kinematic  equations, 

«  =  r  (2.5) 
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1/  =  ain't 


(2.6) 


where  $  is  the  vehicle  heading,  and  y  is  the  cross  track  error  off  a  desired  straight 
line  path. 

B.  COMPENSATOR  DESIGN 

In  control  theory  it  is  known  that  the  eigenvalues  of  the  controller  are  not 
affected  by  the  eigenvalues  of  the  observer.  This  allows  us  to  design  the  controller 
and  observer  separately  which  is  known  as  the  separation  principle.  The  combination 
is  called  a  compensator. 

Equations  (2.3),  (2.5),  and  (2.6)  govern  the  steering  control  of  the  model  used 
in  this  section.  The  control  law  can  be  expressed  as, 

«  =  (2.7) 

where  around  the  nominal  state  '5  =  r  =  y  =  0we  have 

So  =  K*9  +  Krr  +  K„y  (2.8) 

6  is  the  rudder  angle  and  K «,  Kr,  and  Ky  are  the  control  gains  of  the  system. 
The  linear  control  law  is  V  The  rudder  angle  6  is  defined  by  a  hyperbolic  tangent 
function  to  include  the  saturation  to  our  problem  as  shown  in  Figure  2.1.  Saturation 
occurs  at  S^,  which  is  the  saturation  limit  generally  taken  as  0.4  rad. 

The  linearized  form  of  equations  of  motions  in  the  vicinity  of  ^  =  r  =  y  —  0 

are, 


9  =  r 

(2.9) 

r  =  ar  +  b6o 

(2.10) 

y  =  9 

(2.11) 
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Figure  2.1:  Saturation  in  S. 


These  equations  can  be  expressed  in  state  space  form  as 


X  =  AX  4-  Bu 


where 


is  the  state  vector, 


X  = 


r  * 

r 

y  J 


A  = 


0  1  0 
0  a  0 
1  0  0 


is  the  open  loop  dynamics  matrix  and 

B  = 

is  the  control  distribution  vector. 

The  observer  equations  are 


0 

b 

0 


X  =  AX  +  Bu  +  L(Y-  CX) 


(2.12) 


(2.13) 


where  X  is  the  estimate  of  X ,  Y  is  the  output  of  the  system  Y  »  y,  and  C  is  the 
sensor  vector  C  =  [0  0  1]. 


5 


The  error  in  the  estimate  of  X  is  defined  by 


X  =  X  -  X 


Using  equations  (2.12),  (2.13),  and  (2.14)  we  can  obtain 

k  =  (A  -LC)X 


We  can  rewrite  equations  (2.9),  (2.10),  and  (2.11)  in  the  form  of 


‘  0  1  0  ¥ 

=  0  a  0  r 

.  1  0  0  y 


O' 

+  b  6 
0 


The  observer  gains  are, 


After  performing  the  matrix  operations  we  obtain 


Y  =  [0  0  1]  r 

.  y 


r* 

L  —  ir 


*  roi-ftir* 

X=  ?  =  0  a  -lT  f 

y  1  0  y 


Using  equation  (2.13)  we  can  rewrite  equation  (2.8)  as  follows, 
So  =  /?*(¥  -  *)  +  Kr(r  -  f)  +  K,(y  -  y) 


Finally,  we  can  write  our  compensator  equations  in  the  form 


0  1  0  0  0  0 

0  a  0  bcK*  bKr  bK, 

1  0  0  0  0  0 

0  0  0  0  1  -f* 

0  0  0  0  a  -tr 

0  0  0  1  0 


(2.14) 


(2.15) 


(2.16) 


(2.17) 


(2.18) 


(2.19) 
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If  we  look  at  the  matrix  carefully  we  will  see  that  it  is  in  the  form 

’  X  _  A-BK  BK  X  ' 
jC  ~  0  A  -  LC  X 

which  has  the  following  characteristic  equation, 

det[A  -BK- si]  det(A  -  LC  -  si]  =  0 

This  indicates  that  the  dynamics  of  the  observer  are  completely  independent  of  the 
dynamics  (eigenvalues)  of  the  controller.  Thus  K  and  L  can  be  designed  separately. 

C.  CALCULATION  OF  CONTROL  GAINS 

A  is  the  Jacobian  matrix  of  the  system 

’  0  1  O' 

A  —  bK*  a  +  bKr  bKy 
1  0  0 

The  characteristic  equation  of  the  matrix  A  is 

A 3  -  (a  +  bKr)X2  -  bK*  A  -bKy  =  0 
If  the  desired  characteristic  equation  has  the  general  form 

A3  +  orjA3  +  atA  +  ctQ  —  0 
the  control  gains  can  be  found  as 

K  -  ax 

K*  =  "T 

TS  (*2  +  a 

Kr  - - 6“ 

j y  _  °0 

“  b 

The  desired  characteristic  equation  can  be  written  with  respect  to  the  desired  natural 
frequency  and  some  optimum  coefficients.  The  ITAE  criterion  for  a  third  order 
equation  is 

s3  +  1.75t»n«2  +  2.15u>3s  +  =  0 
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where  wn  is  the  desired  controller  natural  frequency. 

Therefore  the  control  gains  can  be  calculated  for  a  given  natural  frequency,  as 

at  =  2.15t0* 


Qj  =  1.75to* 
ao  -  wl 


D.  CALCULATION  OF  OBSERVER  GAINS 


If  we  define  A  as  the  Jacobian  matrix  of  the  system 

0  1  -*#' 

A=  0  a  -tT 
.  1  0  -t,  . 

the  characteristic  equation  of  the  matrix  A  is 


A3  + (*, -a)  +  (** -«*,)A  + (*r -a/*)  =  0 


If  the  desired  characteristic  equation  has  the  general  form 

A3  +  7a*a  +  7i^  +  Tto  *  0 
the  observer  gains  can  be  found  as 

4  =  a  +  7a 
ft  ■  <*4  +  7i 
It  =  a/* +7 o 


Applying  the  ITAE  criteria,  observer  gains  can  be  calculated  for  a  given  natural 
frequency  as 

7,  =  2.15u£0 
72  =  1.75u?„o 

7o  =  u>j|o 
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Normalized  Amp  (  y*  wii~  3  ) 


where  Wno  is  the  observer  natural  frequency. 


E.  CHARACTERISTICS  OF  ITAE  CRITERIA 

In  the  calculations  of  gains  we  applied  the  ITAE  criteria.  If  we  look  at  Figure 
2.2  for  the  step  response  of  ITAE,  we  see  that  the  response  gets  faster  as  the  natural 
frequency  increases.  For  example,  the  settling  time  is  10  normalized  seconds,  or  10 
seconds  for  wn  =  1,  1  second  for  wn  —  10,  and  so  on. 


STEP  RESPONSE  of  ITAE 
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III.  HOPF  BIFURCATION 


A.  INTRODUCTION 

An  important  quantity  in  assessing  the  robustness  of  a  particular  control  law 
design  to  parameter  variations  and  unmodeled  dynamics  is  the  gain  margin.  This  is 
defined  as  the  extent  to  which  changes  can  be  inflicted  on  the  system  gain  without 
loss  in  stability,  to  this  end,  we  assume  that  the  heading  error  gain  K+  is  multiplied 
by  a  constant  C.  By  definition,  of  Hopf  bifurcation  occurs  when  a  pair  of  complex 
conjugate  eigenvalues  cross  into  the  right  hand  half-plane.  When  this  occurs  the 
system  will  deviate  from  a  steady  solution  in  an  oscillatory  manner.  This  deviation 
can  be  either  supercritical  or  subcritical  [Seydel  (1988)].  As  the  parameter  C  crosses 
the  critical  value,  one  pair  of  complex  conjugate  eigenvalues  of  the  linear  system 
matrix  crosses  transversely  the  imaginary  axis.  Locally,  as  C  approaches  C^,  the 
periodic  solutions  are  located  on  the  two  dimensional  Eudedian  plane  spanned  by  the 
eigenvectors  of  the  Jacobian  matrix  of  the  system  which  corresponds  to  the  critical 
pair  of  eigenvalues.  In  this  chap'  -  stability  properties  of  the  periodic  solutions  are 
established.  In  order  to  establish  those  properties  the  main  nonlinear  terms  that 
dominate  the  system  are  isolated.  Center  manifold  theory  is  used  to  reduce  the 
flow  to  a  two  dimensional  manifold.  The  method  of  averaging  is  then  applied  to  the 
reduced  system. 

The  critical  value  of  c  for  stability  of  straight  line  motion  remains  the  same  as 
[Oral  (1993)],  which  is 

C^t  —  0.2658 
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This  is  because  the  dynamics  of  the  controller  are  independent  from  the  dynamics 
of  the  observer  as  explained  in  Chapter  II. 


B.  THIRD  ORDER  EXPANSIONS  OF  THE  SYSTEM  EQUATIONS 
1.  Perturbation  in  A# 

In  the  previous  chapter  we  worked  on  the  linear  system.  Now  we  are  go¬ 
ing  to  introduce  the  nonlinear  terms  to  our  compensator.  In  this  case  the  equations 
of  motion  are 


4f  =  r 

r  =  or  +  aar*  +  W 
y  =  sin® 


where 


=  CK,C»  -  4)  +  K,(r  -  f)  +  K,(y  -  y) 


(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 


or  in  compact  form, 

X  -  f(x) ,  X  =  [#,  r,  y,  *,  f,  »]r  (3.6) 

This  system  can  be  written  in  the  form 

X  =  AX  +  g(x)  (3.7) 

A  is  the  Jacobian  matrix  of  f(x)  evaluated  at  X  =  0,  and  y(x)  contains  all  nonlinear 
terms  of  Equations  (3.1),  (3.2),  and  (3.3).  Taylor  expansion  of  the  nonlinear  terms 
about  the  equilibrium,  where  we  keep  the  first  non-vanishing  nonlinear  coefficients 
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only,  gives 


ain«  = 

(3.8) 

(3.9) 

Substitution  of  Equations  (3.8)  and  (3.9)  into  Equations  (3.1),  (3.2),  and  (3.3)  gives 
us  the  A  matrix  in  Equation  (3.7)  as  follows, 


A  = 


0 

bcK * 
1 
0 
0 
0 


The  nonlinear  parts  are. 


1 

a  +  bKr 
0 
0 
0 
0 


9(x)  = 


0 

bK, 

0 

0 

0 

0 


0 

-bcK* 

0 

0 

0 

1 


0 

-bKr 

0 

1 

a 

0 


1*3 

6 

0 

0 

0 


0 

-bK, 

0 

-fr 


(3.10) 


(3.11) 


If  we  introduce  the  transformation  matrix  ( T )  of  eigenvectors  of  A  evaluated  at  the 
bifurcation  point, 


T  =  KJ  *»j  =  1,2, 3, 4, 5,6  (3.12) 

T-1  =  KJ  *,J  =  1,2, 3, 4, 5, 6  (3.13) 

the  linear  change  of  coordinates, 

x  =  Tz,  z  =  T-'x  (3.14) 

transforms  system  (3.7)  into  its  normal  form 

i  =  T-'ATz  +  T-lg(Tz)  (3.15) 
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With  this  transformation  we  get, 


T~lAT  = 


0  -u*>  0  0  0  0 

Wo  0  0  0  0  0 

0  0  Px  0  0  0 

0  0  0  Pa  0  0 

0  0  0  0  Pa  0 

0  0  0  0  0  PA 


(3.16) 


Using  center  manifold  theory  we  get,  as  shown  in  [Chow  and  Mallet-Paret  (1977)], 

0 


9  = 


^21*1  +  ^22*a  +  ^23*1*2  +  f 24  *1*2 
^31*1  +  ^32*f  +  ^33*1*2  +  ^34*1*3 


(3-17) 


0 
0 
0 

Substitution  of  Equation  (3.14)  into  Equation  (3.7)  yields, 

ii  =  -t0o*2  +  ruzi  +  ruZjZa  +  rx3zxz\  -f  r1Az\  (3.18) 

ia  =  wozi  +  rai^f  +  ^2*23  +  raa2|2a  +  r24z£  (3.19) 

2.  Integral  Averaging 

We  write  Equations  (3.18)  and  (3.19)  in  the  form 


Z\  m  —XVqZ2  +  Fi(Zi,Z2),  (3.20) 

ia  =  woZj  +  Pa(2j,2a)  (3.21) 

If  we  introduce  polar  coordinates  in  the  form, 

Z\  —  RcosO  ,  23  =  Rs'mO  (3.22) 

Equations  (3.20),  (3.21)  result  in 

R  =  Fi(i2, 0)  cos  0  +  Fi(R,  0)  sin  $  (3.23) 

R0  —  wqR  +  J*a(J2, 0 )  c  ’  $  —  Fi(R,  0)  sin  $  (3.24) 
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Equation  (3.23)  yields 


R  =  P(tf)/?3 


(3.25) 


where  P(0)  is  a  2x-periodic  function  in  the  angular  coordinate  0.  If  Equation  (3.25) 
is  averaged  over  one  cycle  in  0,  we  get  an  equation  with  constant  coefficients, 


R  =  KR3 


(3.26) 


where, 

(3.27) 

Equation  (3.27)  is  simplified  after  evaluation  of  the  integral  as, 

K  —  g  {3ru  +  r13  +  rn  +  3r24]  (3.28) 

where  the  coefficients  are  as  follows, 

rii  =  1*13/31  + 1*13/31 
1*13  —  1*13/34  + 1*13/34 
r33  —  1*23/33  + 1*33/33 
1*34  —  1*33/33  + 1*33/33 


where  the  nu,  n13,  n22,  and  1133  are  the  elements  of  the  inverse  of  transformation 
matrix  T.  The  values  of  the  coefficients  Zy  are  given  by  the  following  expressions 

In  =  0317*21  -  bt  [c3A*m4j  +  A?m3 2  +  Zc2 K\KTnx\xmix 

+3C2  KlK^mhm^i  +  ZcK^K^m^m^  +  ZcK^Klm\xmAX  +  ZKrKlrn%xmn 
+3K?Kvmllmti  +  6cA'*A'r/fym41m51m61] 
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in  =  [c3A'|m^  +  ^m|3  + +  3caA’J/irrmJ2m5a 

+Zc7K\  Ktm\2mn  +  ZcK^K}m\2mn  +  ZcK^K^tn^mn  +  3ArrA'^m^3m52 
+3/frJ^rojjni<j  +  6c/ffAfA'yffi4]n)i2ni(]| 

/as  =  3a3m21  6/  j^c3 K\Tn\^Tn42  +  3if*nijjrnj2  +  3ifjnijjW® 

+3c7KlKr  (m^msj  +  2m4im42m81)  +  3 <?K\KV  (m^mej  +  2m4im4jm«i) 

+ZcKtK7  {rn\xm4 a  +  2m5im5jm4iJ  +  3 cK+K7  1m42  +  2m4im4]m4ij 

+ZKrK*  (mJjUljj  +  21719!  7716217151 )  +  ZK*Kt  (mjj  17192  +  2msi  1715217161^ 

+6cKlKrKf  (m4i  17151 17162  +  (l714l77l52  +  7715217151)  77161 ) J 

^24  =  3a3m^n4,  —  6/  |3c® K\m4\rn%2  +  3/iT®»7i5i77i|a  +  ZK^m^rv^ 

+Zc?KlKr  (m^m51  +  2171517715277152]  +  3 <?K\KW  (mjj 77191  +  2771511714317162) 

+ZcK*K?  (771I37715,  +  2771,1 7715217152)  +  ZcKiK7  +  2melmii2m47) 

+SKrKl  (m^TTlsi  +  2 77161 17162 77152)  +  ZK*K9  ^771  j2 17161  +  2l7l5i 771 52 77192) 

+6cKiKrKv  (771521715217161  +  (17151 77152  +  7715217151)  77192)! 

/  —  ^m3 

*31  —  —  gmH 

/  _  1  3 

*32  —  0*7*12 

,  1  3 

*33  =  —  2mnmw 

#  1  3 

*34  =  — -771ii  771 13 

*«  =  3C 
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C.  RESULTS 


The  value  of  K  is  important  for  us  to  determine  whether  the  bifurcation  is 
supercritical  ( K  <  0)  or  subcritical  (K  >  0).  In  this  study  we  wanted  to  see  the 
effects  of  observer  dynamics  to  our  system,  especially  the  value  of  K.  To  do  that 
we  used  the  Fortran  code  (Appendix  A)  for  the  numerical  results.  The  result  was 
that  the  value  of  K  was  not  affected  by  changes  in  the  observer  natural  frequency. 
The  reason  for  this  can  be  traced  back  to  the  definition  of  K ,  Equation  (3.28).  It 
can  be  seen  that  in  the  definitions  for  ri;  and  lij  only  the  first  two  eigenvectors  mtl, 
m,2  for  t  =  1, 2, ...  6  of  the  A  matrix,  Equation  (3.10)  appear.  Therefore,  we  have 
to  show  that  these  remain  the  same  for  all  observer  natural  frequencies. 

This  A  matrix,  Equation  (3.10),  actually  consists  of  4  block  matrices,  each 
3x3,  which  are  the  same  as  shown  in  Equation  (2.22).  Let  us  denote  the  A  matrix 
as  follows, 


A  = 


f  A  B 
0  C 


The  eigenvalues  of  A  can  be  computed  by 


A -XI  B 
0  C-XI 


=  0 


or 

\A  -  AJ|  •  \C  -  A/|  =  0 


We  can  group  the  eigenvalues  in  two  different  groups:  Alt,  for  t  =  1,2,3  are  the 
eigenvalues  of  A  an  A2>,  for  ?  =  1,2,3  are  the  eigenvalues  of  C.  The  eigenvectors 
associated  with  these  eigenvalues  can  be  found  as  follows. 

For  A  =  Alti, 


A  -  A UI 

B 

Vl  _ 

'  0  ' 

0 

C  -  A UI 

Uj  J  " 

0 
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and 

[A  -  +  [B)M  =  0 

MN  +  P-iuflW  -  o 

Since  Aiit  is  an  eigenvalue  of  A  and  the  eigenvalues  of  A  and  C  are  distinct,  the 
matrix  fC  -  Ai,,7]  is  nonsingular  which  means  that  [wj]  =  0.  Therefore,  we  get 

[A  ~  Ax,/1M  =  0 

which  means  that  Vj  is  an  eigenvector  of  A.  Therefore,  the  first  three  eigenvectors  of 
A  are  the  same  as  the  eigenvectors  of  A  and  they  are  independent  of  the  dynamics 
of  the  observer.  Of  course,  the  remaining  three  eigenvectors  of  A  depend  on  the 
oberver  natural  frequency,  but,  as  we  pointed  out  earlier,  none  of  these  appear  in 
the  definition  of  the  nonlinear  stability  coefficient  K. 
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IV.  COMPENSATOR  IN  A  DIFFERENT  BASIS 

A.  CRITICAL  VALUE  OF  C 

If  we  look  at  Equation  (3.6),  we  can  see  that  the  basis  for  our  system  was 
(X,  A].  Now  we  are  going  to  represent  our  system  in  (A,  A')  basis  where  A  is 
the  estimate  of  X.  In  this  compensator  basis  the  critical  value  of  C  in  Equation 
(4.3)  is  no  longer  constant.  Therefore,  we  used  a  Fortran  code  (Appendix  B)  to 
calculate  the  critical  C  values  for  different  observer  natural  frequencies.  A  plot 
of  these  critical  C  values  for  different  observer  natural  frequencies  can  be  seen  in 
Fig.  4.1.  The  observer  natural  frequencies  are  in  nondimensional  seconds  whereas 
the  control  natural  frequencies  are  normalized  with  respect  to  the  corresponding 
observer  natural  frequencies.  The  system  is  unstable  for  values  of  C  less  than  the 
critical  value. 

B.  THIRD  ORDER  EXPANSIONS  OF  THE  SYSTEM  EQUATIONS 

1.  Perturbation  in  A* 

In  the  previous  chapters  we  worked  on  the  [ A,  A]  basis.  Now  we  are 
going  to  represent  our  system  in  the  new  basis,  which  is  (A,  A],  where  A  is  the 
estimate  of  A.  Our  equations  of  motion  were, 

4/  =  r 

r  =  ar  +  a3r3  +  b6  (4.1) 

y  =  sin  ¥ 
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WN(nonnalized) 


Figure  4.1:  Ca *  versus  natural  frequency  for  K 

where 

6  =  ^tanh^j 
So  =  CKi*  +  Krr  +  Kvy 

or  in  compact  form, 

X  */(*)»  A’  =  ['P,r,y,«,r,y]T 
This  system  can  be  written  in  the  form 


A  is  the  Jacobian  matrix  of  f(x)  evaluated  at  X  =  0,  and  $(x)  contains  all  non¬ 
linear  terms  of  Equation  (4.1).  Taylor  expansion  of  the  nonlinear  terms  about  the 
equilibrium,  where  we  keep  the  first  non-vanishing  nonlinear  coefficients  only,  gives 

sin#  =  #-|#3  (4.6) 

*  -  ‘-fib*  <4-7> 

Substitution  of  Equations  (4.6)  and  (4.7)  into  Equation  (4.1)  gives  us  the  A  matrix 

in  Equation  (4.5)  as  follows 

0  1  0  0  0  0 

0  a  0  bcK*  bKr  bK , 

.  1  0  0  0  0  0  /.  q\ 

A=  o  0  /*  0  1  -/#  (4‘8) 

0  0  ir  bcK *  bKr  -It  +  bK, 


0  0  L 


The  nonlinear  parts  are, 


9(x)  = 


i 

_±<p3 


b  6> 
~3&T° 


If  we  introduce  the  transformation  matrix  (T)  of  eigenvectors  of  A  evaluated  at  the 
bifurcation  point, 


T  =  [m„]  i,;  =  1,2, 3, 4, 5, 6 

T-1  =  K]  i,j  =  1,2, 3,4, 5, 6 


(4.10) 

(4.11) 


the  linear  change  of  coordinates, 


x  —  Tz  ,  z  —  T~l: 


(4.12) 
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transforms  system  (4.5)  into  its  normal  form 


i  =  Tx  AT z  +  T~lg(Tz) 


(4.13) 


At  the  bifurcation  point 


t~xat  = 


0  -too  0  0  0  0 

too  0  0  0  0  0 

0  0  Px  0  0  0 

0  0  0  Pi  0  0 

0  0  0  0  P9  0 

0  0  0  0  0  PA 


(4.14) 


with  too  >  0  and  P,  <  0. 


Using  center  manifold  theory  we  get,  as  shown  in  [Chow  and  Mallet-Paret 


(1977)], 


*21*1  +  fa*?  +  *23*1*2  +  *24*1*2 
*31*1  +  *32*2  +  *33*1*2  +  *34*1*2 
9  0 

*51*1  +  *52*f  +  *53*1*2  +  *54*1*2 
0 

Substitution  of  Equation  (4.12)  into  Equation  (4.5)  yields, 

*1  =  -mo*2  +  rnzf  +  ri3S?za  +  r^zizl  +  rXAz\ 
*2  —  *00*1  +  r21*l  +  **22*1*2  +  *"23*1*2  +  r24*2 


(4.15) 


(4.16) 

(4.17) 


where  the  terms  nj  are  evaluated  by  a  Fortran  code  (Appendix  C). 

For  values  of  C  close  to  its  critical  value.  Equations  (4.16)  and  (4.17) 

become, 

i\  =  a!zzx  -  (too  +  *o'c)*2  +  ruz\  +  ri2*?*2  +  risZisf  +  ruz%  (4.18) 
*2  =  (wo  +  to'e)zj  +  oltzi  +  rnz\  +  r^z^zj  +  r23*i*|  +  **34*f  (4.19) 

where  e  is  the  difference  between  C  and  its  critical  value  C*.  The  terms  a'  and  to* 
denote  the  derivative  of  the  real  and  imaginary  part,  respectively,  of  the  critical  pair 
of  the  eigenvalues  with  respect  to  C  evaluated  at  C* 
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2.  Integral  Averaging 

We  write  Equations  (4.18)  and  (4.19)  in  the  form, 

it  =  a'ezi  -  (tuo  +  tn'e)z3  +  /\(zt,  z3)  (4.20) 

z3  =  (t£H>  +  to'c)z i  +  oltz-i  +  F3(zu  z3)  (4.21) 

If  we  introduce  polar  coordinates  in  the  form, 

Zi  =  Rax  $  ,  Zi  =  Ram  6  (4.22) 

Equations  (4.20),  (4.21)  result  in 

R  =  a'eR  +  Fl(R>O)co60  +  F3(R,0)  am  0  (4.23) 

R0  =  (uto  +  u/e)/2  +  .F3(/2,0)co8  0  —  Ft(R,0)  sin0  (4.24) 

Equation  (4.23)  yields 

R  =  a'eR  +  P(0)R?  (4.25) 


where  P(0)  is  a  2x-periodic  function  in  the  angular  coordinate  0.  If  Equation  (4.25) 
is  averaged  over  one  cycle  in  0  [Chow  and  Mallet-Paret  (1977)],  we  get  an  equation 
with  constant  coefficients, 

R  =  a'eR+KR3  (4.26) 

where 

K  =  h[mie  (4-27) 

Equation  (4.27)  is  simplified  after  evaluation  of  the  integral  as, 

K  —  g  [3rtt  ■+■  rj3  +  +  3tm]  (4.28) 

where  the  coefficients  are  as  follows 
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rll  =  n12^21  +  n  13^31  +  n15^51 


T|3  =  ”12^24  +  »13^34  +  n15^M 

rn  —  1133/33  +  7*23/33  +  7135/53 
rj4  =  n  23/22  +  7*23/32  +  n  25/52 

where  the  nja,  nj3,  nj5, 7*22,  1*23,  and  nas  are  the  elements  of  the  inverse  of  trans¬ 
formation  matrix  T.  The  values  of  the  coefficients  /21,  /3a,  /as*  Z24,  /31,  /aa,  /3s,  and 
/s4  are  the  same  as  in  Chapter  III.  The  values  of  the  other  /,j  coefficients  are  given 
by  the  following  expressions: 

/si  =  -6,  [c3/^  +  Krmh  +  +  &KlKTm\xmtl 

+3<?KlK,mllm1il  +  3cA'*A'r2m^1m41  +  3cK*Klrr&mAX  +  3/Tr  ffjn&m,! 
+3K*Ktm\lmti  +  6cKtKrKym4XmsXm^x^ 

la  =  -64  [c3 A'Jm’j  +  K*n&  +  Kjrn*,  + 

+3c*  fif *  Kvm\2Tna  +  3ctf*/irram5a"*43  +  ZcKiKlrr&m*  +  3KrK}n&mn 

+ZK*KyTn\2m€2  +  6cA’*/irrA’,m43m53m<{3j 

/*3  =  —be  |3c? K%m\lm42  +  ZKfm\xm$2  +  ZK^rn^im^ 

+3 <?K\Kr  (m^msa  +  2m41m4amSi )  +  Z<?K\Ky  (m41m«a  +  2m41m4am«i) 
+3cKiK?  (m|1m42  +  2m51m52m41)  +  ZcK^K]  (mj,m42  +  2m«1m«am4i) 
+3 KrK\  (mj,m| 3  +  +  ZK*Ky  (m,tm«2  +  2m51m82m61) 

+6 cK*KrKt  (m41m5im«2  +  (m41m5a  +  m4amsi)  171*1)] 
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/m  —  —ftr  [3c* +  3/f*msimJ3  +  ZK^m^rn^ 

+3c* K\Kr  (jn^mn  +  +  3c3 K\K9  ^m^mn  +  2n>4i TO42 m<jj 

+3 cKiK*  (m33m4i  +  +  ZcK9K*  (rn^m^  +  2m«me2m4]j 

+ZKrKl  (mJa m51  +  2m61me3m53)  +  ZK?Kt  (m^3m«  +  2m51m53m63) 
+ZcK+KrKv  (m^msjm*!  +  (m4ini|3  +  m«)| 

C.  RESULTS 

Existence  and  stability  of  limit  cycles  can  be  determined  by  analyzing  the 
equilibrium  points  of  the  averaged  Equation  (4.26),  which  correspond  to  periodic 
solutions  in  Zj,  z3  as  can  be  seen  from  Equation  (4.22).  FYom  Equation  (4.26)  we 
can  easily  see  that: 

1.  If  a'  >  0,  then 

(a)  if  K  >  0,  then  unstable  periodic  solutions  co-exist  with  the  stable  equi¬ 
librium  for  e  <  0,  and 

(b)  if  K  <  0,  then  stable  periodic  solutions  co-exist  with  the  unstable  equi¬ 
librium  for  e  >  0. 

2.  If  ct  <  0,  then 

(a)  if  K  >  0,  then  unstable  periodic  solutions  co-exist  with  the  stable  equi¬ 
librium  for  e  >  0,  and 

(b)  if  K  <  0,  then  stable  periodic  solutions  co-exist  with  the  unstable  equi¬ 
librium  for  e  <  0. 
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We  refer  to  K  <  0  as  the  supercritical,  and  K  >  0  as  the  subcritical  PAH  bi¬ 
furcation.  In  the  supercritical  case,  after  the  equilibrium  state  loses  its  stability 
the  system  converges  to  a  stable  periodic  solution  with  amplitude  which  increases 
continuously  as  the  difference  e  is  increased. 

In  the  subcritical  case,  however,  before  the  equilibrium  state  loses  stability,  its 
domain  of  attraction  becomes  very  small  since  it  is  bounded  by  the  amplitudes  of  the 
unstable  limit  cycles.  In  such  a  case,  an  initial  disturbance  of  sufficient  magnitude 
can  throw  the  system  off  the  nominal  path  even  before  its  domain  of  attraction 
has  completely  shrunk  to  zero.  As  the  nominal  equilibrium  becomes  unstable,  the 
system  jumps  to  a  different  state  of  motion  with  a  locally,  at  e  —  0,  discontinuous 
increase  in  the  amplitude  [Papoulias  (1993)]. 

In  our  case,  the  value  of  ol  is  always  negative,  which  can  be  seen  easily  from 
Figure  4.1.  If  we  look  at  the  nature  of  the  curve  for  the  critical  value  of  C  for 
different  natural  frequencies,  we  will  see  that  as  the  value  of  critical  C  decreases  for 
the  same  natural  frequency,  the  system  becomes  unstable. 

After  using  a  Fortran  code  (Appendix  C)  we  observed  that  the  nonlinear  sta¬ 
bility  coefficient  K  depends  on  observer  dynamics.  Figure  4.2  shows  that  for  a  given 
control  design,  the  observer  must  be  as  responsive  as  possible  to  ensure  negative  K 
(stable  limit  cycle).  On  the  other  hand,  for  a  given  observer  design,  if  the  control 
law  is  too  slow  we  get  subcritical  behavior  (unstable  limit  cycle). 
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Figure  4.2:  Kk+  versus  wn  for  different  observer  wn. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

An  investigation  of  the  nonlinear  dynamic  response  characteristics  of  a  marine 
vehicle  has  been  presented.  Particular  emphasis  in  this  work  was  placed  on  analyzing 
the  effects  of  observer  design  on  system  response  after  initial  loss  of  stability  of 
straight  line  motion.  Bifurcation  theory  techniques  were  utilized  in  order  to  assess 
that  behavior.  The  main  conclusions  of  this  work  can  be  summarized  as  follows. 

1.  There  exists  a  critical  point  for  a  certain  combination  of  system  gains  and 
system  parameters  for  stability  of  straight  line  motion.  The  loss  of  stability 
occurs  generically  in  the  form  of  Poincare-Andronov-Hopf  bifurcations.  As 
the  parameter  crosses  its  critical  value,  a  family  of  periodic  orbits,  self  sus¬ 
tained  oscillations  develops.  Center  manifold  reduction  and  integral  averaging 
techniques  were  used  in  order  to  establish  the  direction  of  the  bifurcation  and 
stability  of  the  resulting  periodic  solutions  [Papoulias,  Oral  (1993)]. 

2.  For  [X,  X]  basis  the  critical  point  does  not  depend  on  the  observer  dynamics 
(separation  principle).  The  nonlinear  stability  coefficient  K  was  not  influenced 
by  observer  dynamics.  The  previous  reduction  process  shows  that  K  depends 
on  the  first  two  eigenvectors  of  the  6x6  matrix  A.  Matrix  algebra  shows  that 
these  eigenvectors  are  associated  only  with  the  controller  dynamics. 

A 

3.  For  [X,  X]  basis  the  critical  value  depends  on  observer  dynamics.  For  a  given 
control  design,  the  observer  must  be  as  responsive  as  possible  to  maximize  the 
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region  of  stability.  On  the  other  hand,  for  a  given  observer  design,  the  control 
must  be  as  slow  as  possible  to  maximize  region  of  stability. 

4.  The  nonlinear  stability  coefficient  K  depends  on  observer  dynamics  for  this 
basis.  For  a  given  control  design,  the  observer  must  be  as  responsive  as  possible 
to  ensure  negative  K  (stable  limit  cycles).  In  this  benign  loss  of  stability 
the  resulting  periodic  solutions  are  continuous  single-valued  functions  of  the 
parameter  distance  from  its  critical  value.  On  the  other  hand,  for  a  given 
observer  design,  if  the  control  law  is  too  slow  we  get  subcritical  behavior 
(unstable  limit  cycles).  In  such  a  case,  the  periodic  solutions  develop  with 
what  appears  to  be  a  discontinuous  increase  in  the  amplitude  of  oscillations 
[Papoulias,  Oral  (1993)]. 

B.  RECOMMENDATIONS 

The  differences  between  the  two  bases  with  respect  to  robustness  properties  of 
the  system  have  to  be  analyzed. 
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APPENDIX  A 

HOPF  BIFURCATION  PROGRAM  FOR  [X,X]  BASIS 


PROGRAM  HTKPSI 
HOPF  BIFURCATIONS 
NOMOTO'S  FIRST  ORDER  MODEL 

CALCULATIONS  FOR  K  AND  CCRITICAL  IF  KPSI  CHANGES  W/  C 
C234567 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C 

DOUBLE  PRECISION  K1 , K2 , K, LPSI , LY, LR,  IZ,L, 

&  MA3S,NV,NRiNVDOT,NRDOn?,NDRS,NDRB,  KPSI,KR,KY,K3, 

&  L21 , L22 , L23 , L24 , L3 1 , L32 , L3  3 , L34 , L51 , L52 , L53 ,  L54 , 

&  Ml  1,M12,M13/M14jM15,M1S/M21,M22,M23,M24,M25,M26, 

&  M31 ,  M32 ,  M33 ,  M34 ,  M3  5 ,  M3  6 ,  M41 ,  M42 ,  M43 ,  M44 ,  M45 ,  M46 , 

&  M51 ,  M52  ,  M53 ,  M54 ,  M55 ,  M56 ,  M61 ,  M62 ,  M63 ,  M64 ,  M65 ,  M66 , 

&  Nil,Nl2,N13,N14,N15,N16,N21,N22,N23,N24,N25,N26, 

&  N31,N32,N33,N34,N35,N36,N41,N42  ,N43,N44,N45,N46, 

&  N51 ,  N52 ,  N53 ,  N54 ,  N55 ,  N56 ,  N61 ,  N62 ,  N63 ,  N64 ,  N65 ,  N66 

C 

DIMENSION  AMAT  (6,6)  ,T(6,6)  ,  TINV(  6, 6)  ,  FV1  ( 6 )  ,  IV1  { 6 )  ,YYY(6,6) 
DIMENSION  WR(  6)  ,WI(6>  ,TSAVE(6,6)  ,TLUD(6,6)  ,  IVLUD{6)  ,  SVLUD{6) 
DIMENSION  ASAVE  (6,6)  ,Al(6f6)  ,  A2  ( 6 , 6 ) 

OPEN  (11, FILE= ' AKPSI . MAT ' , STATUS = ' NEW ' ) 

C 

WEIGHT=435 . 0 

IZ  =45.0 

L  =7.3 

RHO  =1 . 94 

G  =32.2 

XG  =0.0104 

MASS  =WEIGHT/G 

MASS  =MASS/(0.5*RHO*L**3) 

IZ  =IZ/ (0.5*RHO*L**5) 

XG  =XG/L 
YRDOT  =-0.00000 
YVDOT  =-0.03430 
YR  =+0.00000 
YV  =-0.10700 
YDRS  =+0.01241 
YDRB  =+0.01241 
NRDOT  =-0.00047 
NVDOT  =-0.00000 
NR  =-0.00390 
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NV  =-0.00000 
NDRS  =-0.337*YDRS 
NDRB  =+0.283*YDRS 
DH  = ( IZ-NRDOT) * (MASS-YVDOT) - 
&  (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 

All= ( { IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DH 
A12=( (IZ-NRDOT) * (-MASS+YR) - 
&  (MASS*XG-YRDOT) * ( -MASS*XG+NR) ) /DH 

A21= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DH 
A22=(  (MASS-YVDOT)  *  ( -MASS*XG+NR)  - 
&  (MASS*XG-NVDOT) * ( -MASS+YR) ) /DH 

Bll= ( (IZ-NRDOT) *YDRS- (MASS*XG-YRDOT) *NDRS ' /DH 
B12= ( (IZ-NRDOT) *YDRB- (MASS*XG-YRDOT) *NDRB) /DH 
B21= { (MASS-YVDOT) *NDRS- (MASS*XG-NVDOT) *YDRS) /DH 
B22= ( (MASS-YVDOT) *NDRB- (MASS*XG-NVDOT) *YDRB) /DH 
B1  =B11-B12 
B2  =B21-B22 
C 

200  WRITE  (*,1004) 

READ  (*,*)  WNMIN , WNMAX , IWN 

INCR=IWN 

WRITE  (*,*)  'ENTER  OBSERVER  WN' 

READ  (*,*)  WNO 
205  WRITE ( * , 1007 ) 

READ  (*,*)  A3 
C  DO  is  Dsat 

50  WRITE  (*,1006) 

READ  ( * ,  * )  DO 

WRITE  (*,1008) 

READ  ( * , * )  CCRIT 

C1=(A11*A22-A21*A12) * (A21*B1-A11*B2 ) 

C2= (A11+A22 ) * (A21*B1-A11*B2 ) +B2* (A11*A22-A21*A12 ) 
C3=-(A21*B1-A11*B2) **2 
A=C1/C2 
B=C3/C2 
C 

DO  1  11=1, INCR 
C 
C 

WN  =WNMIN+ (WNMAX -WNMIN) * (II— 1) / (INCR-1) 

C  print  *,wn 

ALPHA0  =WN*  *  3 
ALPHA1=2 . 15*WN**2 
ALPHA2=1.75*WN 
KPSI=-ALPHA1/B 
KY  =-ALPHA0/B 
KR  =- (ALPHA2+A) /B 
GAMAO  =WNO*  *  3 
GAMA1 =2 . 1 5 *WNO* * 2 
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GAMA2=1 . 75*WNO 
LY=A+GAMA2 
LPSI=A*LY+GAMA1 
LR=A*LPSI+GAMAO 

C234567 

C  A  MATRIX 

AMAT{1, 1)  =0 . 0 

AMAT ( 1 , 2 ) =1 . 0 

AMAT{1 , 3 ) =0 . 0 

AMAT ( 1, 4 ) =0 . 0 

AMAT ( 1 , 5 ) =0 . 0 

AMAT(1, 6) =0 . 0 

AMAT  (2,1) =B*CCRIT*KPSI 

AMAT (2 , 2 ) =A+ (B*KR) 

AMAT (2,3) =B*KY 
AMAT (2,4) =-B*CCRIT*KPSI 
AMAT (2, 5) =-B*KR 
AMAT (2, 6) =-B*KY 
AMAT ( 3 , 1 ) =1 . 0 
AMAT (3 , 2 ) =0 . 0 
AMATO  ,  3 )  =0 . 0 
AMAT(3, 4)=0.0 
AMATO,  5)  =0.0 
AMAT ( 3 , 6 ) =0 . 0 
AMAT (4, 1 ) =0 . 0 
AMAT (4 , 2 ) =0 . 0 
AMAT (4, 3 )  =0 . 0 
AMAT (4, 4) =0 . 0 
AMAT ( 4 , 5 ) =1 . 0 
AMAT (4,6) =-LPSI 
AMAT (5, 1)  =0 . 0 
AMAT(5, 2 ) =0 . 0 
AMAT(5, 3  )  =0 . 0 
AMAT(5, 4) =0 . 0 
AMAT (5,5) =A 
AMAT  (5,6)  =-LR 
AMAT (6, 1) =0 . 0 
AMAT ( 6 , 2 ) =0 . 0 
AMAT (6, 3) =0 .0 
AMAT ( 6, 4)  =1 . 0 
AMATO,  5)  =0 . 0 
AMAT  (6,6)  =-LY 
DO  11  1=1,6 
DO  12  J=l, 6 

ASAVE ( I , J) =AMAT ( I , J) 

12  CONTINUE 

11  CONTINUE 

CALL  RG (6, 6, AMAT,WR, WI, 1, YYY, IV1, FV1, IERR) 
CALL  DSOMEG(IEV,WR,WI, OMEGA, CHECK) 

C  WRITE  (*,*)  IEV 
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WRITE  (60,*)  (WR(IREAL) ,IREAL=1, 6) 

OMEGAO  =OMEGA 
DO  5  1=1,6 

T(I, 1) =YYY (I, IEV) 

T(I, 2) =-YYY (I, IEV+1) 

5  CONTINUE 

IF(IEV.EQ.l)  GO  TO  13 
IF ( IEV. EQ . 2 )  GO  TO  14 
IF (IEV.EQ . 3 )  GO  TO  15 
IF ( IEV. EQ . 4 )  GO  TO  16 
IF (IEV.EQ. 5)  GO  TO  17 
STOP  3004 

13  DO  21  1=1,6 

T(I, 3) =YYY (1,3) 

T ( I , 4 ) =YYY (1,4) 

T ( I , 5 ) =YYY (1,5) 

T(I, 6) =YYY (1,6) 

21  CONTINUE 
GO  TO  30 

14  DO  22  1=1,6 

T ( I , 3 ) =YYY (1,1) 

T  ( 1 , 4 ) =YYY (1,4) 

T ( I , 5 ) =YYY (1,5) 

T ( I , 6 ) =YYY (1,6) 

22  CONTINUE 
GO  TO  30 

15  DO  23  1=1,6 

T  ( 1 , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T ( I , 5 ) =YYY (1,5) 

T ( I , 6 ) =YYY (1,6) 

23  CONTINUE 
GO  TO  30 

16  DO  24  1=1,6 

T ( I , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T ( I , 5 ) =YYY (1,3) 

T ( I , 6 ) =YYY (1,6) 

24  CONTINUE 
GO  TO  30 

17  DO  25  1=1,6 

T(I, 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T ( I , 5 ) =YYY (1,3) 

T(I, 6) =YYY (1,4) 

25  CONTINUE 

30  CONTINUE 

C 

C  NORMALIZATION  OF  THE  CRITICAL  EIGENVECTOR 

C 
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CALL  NORMAL (T) 

INVERT  TRANSFORMATION  MATRIX 

DO  2  1=1,6 
DO  3  J=l, 6 

TINV { I , J) =0 . 0 
TSAVE ( I , J ) =T ( I , J ) 

CONTINUE 

CONTINUE 

CALL  DLUD (6,6, TSAVE, 6 , TLUD, IVLUD) 
DO  4  1=1,6 

IF  ( IVLUD { I )  . EQ . 0 )  STOP  3003 
CONTINUE 

CALL  DILU (6,6, TLUD , IVLUD , SVLUD ) 

DO  8  1=1,6 
DO  9  J=l, 6 

TINV ( I , J ) =TLUD ( I , J ) 

CONTINUE 

CONTINUE 

CHECK  Inv(T)*A*T 


IMULT=1 

IF  (IMULT.EQ.l)  CALL  MULT (TINV, ASAVE,T,A2) 
IF  (IMULT.EQ.0)  STOP  3007 
P1=A2 (1,1) 

P2=A2 (2,2) 

P=A2  (3,3) 

Q=A2 (4,4) 

R=A2 (5,5) 

S=A2 (6,6) 

WRITE  (21,  *)P1,P2,P,Q,R,S 

DEFINITION  OF  Nij 

N11=TINV(1, 1) 

N21=TINV(2 , 1 ) 

N31=TINV (3,1) 

N41=TINV (4,1) 

N51=TINV(5, 1) 

N61=TINV (6,1) 

N12=TINV (1,2) 

N22=TINV(2 , 2 ) 

N32=TINV (3,2) 

N42=TINV (4,2) 

N52=TINV (5,2) 

N62=TINV(6, 2) 

N13=TINV (1,3) 

N23=TINV (2,3) 
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N33=TINV(3,3) 
N43=TINV(4,3) 
N53=TINV (5,3) 
N63=TINV<6, 3) 
N14=TINV(1,4) 
N24=TINV(2,4) 
N34=TINV (3,4) 
N44=TINV(4,4) 
N54=TINV(5, 4) 
N64=TINV{6, 4) 
N15=TINV(1,5) 
N25=TINV(2 , 5) 
N35=TINV (3,5) 
N45=TINV(4, 5) 
N55=TINV(5, 5) 
N65=TINV (6,5) 
N16=TINV (1,6) 
N26=TINV<2,6) 
N36=TINV (3,6) 
N46=TINV(4,6) 
N56=TINV (5,6) 
N66=TINV(6, 6) 

C 

C  DEFINITION  OF  Mij 

C 

M11=T(1, 1) 

M21=T (2 , 1 ) 

M31=T (3 , 1 ) 
M41=T(4, 1) 
M51=T(5, 1) 
M61=T(6, 1) 
M12=T(1, 2 ) 

M22=T(2 , 2 ) 

M32=T(3 , 2 ) 
M42=T(4, 2 ) 
M52=T(5, 2 ) 
M62=T(6, 2 ) 
M13=T(1, 3 ) 
M23=T(2, 3) 

M33=T (3 , 3 ) 
M43=T(4, 3 ) 
M53=T(5, 3 ) 
M63=T(6, 3 ) 
M14=T(1, 4) 

M24=T{2 , 4) 

M34=T ( 3 , 4 ) 
M44=T(4, 4) 
M54=T(5, 4) 
M64=T{6, 4) 
M15=T(1, 5) 
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M25=T(2, 5) 

M35=T(3, 5) 

M45=T(4, 5) 

M55=T(5, 5) 

M65=T(6, 5) 

M16=T (1,6) 

M26=T(2 , 6) 

M36=T ( 3 , 6 ) 

M46=T(4, 6} 

M56=T(5, 6) 

M66=T(6, 6) 

WRITE (70, *)  N11,N12,N13 
WRITE (71,*)  N14,N15,N16 
WRITE (72 , * )  N21,N22,N23 
WRITE (73 , * )  N24,N25,N26 
WRITE (74 , * )  N3 1 , N3 2 , N3 3 
WRITE (75, *)  N34,N35,N36 
WRITE (76, *)  N41,N42,N43 
WRITE (77,*)  N44,N45,N46 
WRITE (78,*)  N51,N52,N53 
WRITE (79 , * )  N54,N55,N56 
WRITE (80, *)  N61,N62,N63 
WRITE (81, *)  N64,N65,N66 
Kl=l . /8 . * ( (ALPHA2**3 ) +ALPHA0 ) / (ALPHA2 ) 

K2=3.*A3-.5* ( ALPHA2 *  *  2 ) / ALPHAO 

K3=l./( (B**2)*(D0**2) ) * (ALPHA2+A) * ( (ALPHA0/ALPHA2 ) + (A**2 ) ) 
K=K1* (K2+K3 ) 

print  *,  wn,k,ccrit 

BL=B/ (3*D0**2) 

C2 3 45 67890123456789012345678 901234567 8901234567890123456789012345 
6789012 

L21=A3*M21**3-BL* (CCRIT**3*KPSI**3*M41**3+KR**3*M51**3+ 

&  KY**3*M61**3+ 

&  3*CCRIT**2*KPSI**2*KR*M41**2*M51+ 

&  3*CCRIT**2*KPSI**2*KY*M41**2*M61+ 

&  3*CCRIT*KPSI*KR**2*M51**2*M41+3*CCRIT*KPSI*KY**2*M61**2*M41+ 
&  3*KR*KY**2*M61**2*M51+3*KR**2*KY*M51**2*M61+ 

&  6*CCRIT*KPSI*KR*KY*M41*M51*M61) 

L22=A3*M22**3-BL* (CCRIT**3*KPSI**3*M42**3+KR**3*M52**3+ 

Sc  KY**3*M62**3+ 

Sc  3  *CCRIT*  *2  *KPSI**2  *KR*M42  *  *2  *M52+ 

&  3  *CCRIT* *  2  *KPSI *  *2  *KY*M42  *  *2  *M62+ 

Sc  3*CCRIT*KPSI*KR**2*M52**2*M42+3*CCRIT*KPSI*KY**2*M62**2*M42+ 
&  3*KR*KY**2*M62**2*M52+3*KR**2*KY*M52**2*M62+ 

Sc  6*CCRIT*KPSI*KR*KY*M42*M52*M62) 

L23=3*A3*M21**2*M22-BL* (3*CCRIT**3*KPSI**3*M41**2*M42+ 

&  3*KR**3*M51**2*M52+ 
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&  3*KY**3*M61**2*M62+ 

&  3  *CCRIT**2*KPSI**2*KR* (M41 *  *  2  *M52+2  *M41 *M42  *M51 )  + 

&  3*CCRIT**2*KPSI**2*KY* (M41**2*M62+2*M41*M42*M61) + 

&  3  *CCRIT* KPSI * KR*  *2  *  < MS 1 *  * 2  *M42 +2  *M51 *M52  *M41 )  + 

&  3*CCRIT*KPSI*KY**2MM61**2*M42+2*M61*M62*M41)+ 

&  3*KR*KY**2* (M61**2*M52+2*M61*M62*M51) + 

&  3*KR**2*KY* (M51**2*M62+2*M51*M52*M61) + 

&  6*CCRIT*KPSI**KR*KYMM41*M51*M62+(M41*M52+M42*M51)  *M61)  ) 

L24=3*A3*M21*M22**2-BL* (3*CCRIT**3*KPSI**3*M41*M42**2+ 

&  3*KR**3*M51*M52**2+ 

&  3*KY**3*M61*M62**2+ 

&  3  *CCRIT**2*KPSI *  *2  *KR* (M42  *  *2  *M51+2  *M41*M42  *M52 )  + 

&  3  *CCRIT* *2  *KPSI*  *2  *KY* (M42  *  *2 *M61+2 *M41  *M42  *M62 )  + 

&  3  *CCRIT*KPSI*KR*  *2  * (M52**2*M41+2*M51*M52*M42 )  + 

&  3  *CCRIT*KPSI *KY*  *2  * (M62  *  *2 *M41+2 *M61 *M62  *M42 )  + 

&  3*KR*KY**2*(M62**2*M51+2*M61*M62*M52)+ 

&  3*KR**2*KY*(M52**2*M61+2*M51*M52*M62)+ 

&  6*CCRIT*KPSI*KR*KY* <M42*M52*M61+ (M41*M52+M42*M51) *M62) ) 

L31=(-1./6.)*M11**3 
L32= (— 1./6. ) *M12**3 
L33= (-1 . /2 . ) *M11**2*M12 
L34= ( -1 . /2 . ) *M11*M12**2 
R11=N12*L21+N13*L31 
R12=N12*L23+N13*L33 
R13=N12*L24+N13*L34 
R14=N12*L22+N13*L32 
R21=N22*L21+N23*L31 
R22=N22*L23+N23*L33 
R23=N22*L24+N23*L34 
R24=N22*L22+N23*L32 
C 

K= (3*R11+R13+R22+3*R24) /8 
WRITE  (11,2001)  WN, K, CCRIT 
1  CONTINUE 
STOP 

1004  FORMAT  ('  ENTER  MIN, MAX,  AND  INCREMENTS  OF 

WN(stepstodivide) ' ) 

1006  FORMAT  ('  ENTER  DELTASAT.') 

1007  FORMAT  ('  ENTER  A3 ' ) 

1008  FORMAT  ('  ENTER  CCRIT') 

2001  FORMAT  (3E15.5) 

END 

SUBROUTINE  DSOMEG ( IJK, WR,WI, OMEGA, CHECK) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  WR ( 6 ) , WI ( 6 ) 

CHECK=-1 . OD+25 
DO  1  1=1,6 

IF  (WR(I) .LT. CHECK)  GO  TO  1 
CHECK=WR ( I ) 
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IJ=I 

1  CONTINUE 

OMEGA = DABS (WI { IJ) ) 

IF  (WI (IJ) .GT.O.DO)  IJK=IJ 


IF  (WIlIJ) .LT.O.DO)  IJK=IJ-1 

RETURN 

END 

SUBROUTINE  NORMAL (T) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  T { 6 , 6 ) , TNOR (6,6) 

CFAC=T(1,1)**2+T(1,2)**2 

IF  (DABS (CFAC) .LE. (l.D-10) )  STOP  4001 

TNOR ( 1 , 1 ) si . DO 

TNOR(2,l)=(T(l,l) *T(2,1)+T(2,2)*T(1,2) ) /CFAC 
TNOR (3 , 1 )  =  (Td,  1)  *T(3 , 1)  +T{3 , 2 )  *T(1, 2 )  ) /CFAC 
TNOR (4,1)=(T(1,1)*T(4,1) +T ( 4 , 2 ) *T (1 , 2 ) ) /CFAC 
TNOR (5,1)=(T(1,1)*T(5,1) +T(5, 2 ) *T(1, 2 ) ) /CFAC 
TNOR(6,l)=(T(l,l)*T(6,l)+T(6,2)*T(l,2) ) /CFAC 
TNOR(l, 2 ) =0 .DO 

TNOR (2,2)=(T{2,2)*T(1/1)“T(2,1)*T(1,2)) /CFAC 
TNOR (3,2)=(T(3,2) *T (1,1) -T ( 3 1 1 ) *T ( 1 , 2 ) ) /CFAC 
TNOR(4,2)=(T(4/2)*T(l,l)-T(4,l)*T(l,2) ) /CFAC 
TNOR (5, 2) * (T(5, 2) *T(1, 1) -T(5, 1) *T(1, 2) ) /CFAC 
TNOR (6,2)=(T(6,2)*T(1,1)“T(6,1)*T(1,2)) /CFAC 
DO  1  1*1,6 
DO  2  J=l, 2 

T(I, J) =TNOR(I, J) 

2  CONTINUE 

1  CONTINUE 
RETURN 
END 

C======*====*===============**===*=============================*= 

SUBROUTINE  MULT (TINV, A, T, A2 ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  TINV (6,6) ,A(6,6),T(6,6),A1(6/6),A2(6,6) 

DO  1  1*1,6 
DO  2  J=l, 6 
A1(I, J)=0.D0 
A2(I,J)=0.D0 

2  CONTINUE 
1  CONTINUE 

DO  3  1=1,6 
DO  4  J=l, 6 
DO  5  K=l, 6 

AKI,  J)=A(I,K)*T(K,J)+A1(I,J) 
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5  CONTINUE 

4  CONTINUE 
3  CONTINUE 

DO  6  1=1,6 
DO  7  J=l, 6 
DO  8  K=l, 6 

A2  (I,  J)  =TINV (I, K)  *A1  (K,  J)  +A2  (I ,  J) 
8  CONTINUE 


7  CONTINUE 
6  CONTINUE 

DO  11  1=1,6 

C  WRITE  (*,101)  (A(I,J),J=1,6) 


11 

CONTINUE 

DO  12  1=1,6 

c 

12 

WRITE  (*,101) 
CONTINUE 

DO  10  1=1,6 

(T(X,J) , J=l, 6) 

C 

10 

WRITE  (*,101) 
CONTINUE 

(A2  (I,  J)  ,  J=l,  6) 

C 

WRITE  (*,101) 
RETURN 

A2(l,l) 

101 

FORMAT  (4E15 
END 

.5) 
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APPENDIX  B 

CRITICAL  VALUE  OF  C  FOR  [X,X]  BASIS 


PROGRAM  NCCRIT 
HOPF  BIFURCATIONS 
NOMOTO'S  FIRST  ORDER  MODEL 

CALCULATIONS  FOR  CCRITICAL 

C234567 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1 , K2 , K, LPSI , LY, LR, IZ, L, 

&  MASS , NV, NR, NVDOT, NRDOT, NDRS , NDRB, KPSI , KR, KY, K3 

C 

DIMENSION  AMAT (6,6) , FV1 ( 6 ) , IV1 ( 6 ) ,YYY(6,6) 

DIMENSION  WR  (6)  ,WI(6)  ,  TSAVE  (6,  6)  ,  TLUD  ( 6,  6)  ,IVLUD(6)  ,SVLUD(6) 
DIMENSION  ASAVE (6,6)  ,A1(6,6)  ,A2  (6,6) 

OPEN  ( 11, FILE* 'CVWN1. RES' , STATUS* ' NEW' ) 

OPEN  (12,FILE='CVWN2 .RES' , STATUS* 'NEW' ) 

OPEN  ( 13 , FILE= ' CVWN3 . RES ' , STATUS* ' NEW ' ) 

WRITE  (*,*)  'ENTER  MIN, MAX,  AND  INCREMENTS  IN  CCRIT' 

READ  (*,*)  CMIN, CMAX, IC 

WRITE  (*,*)  'ENTER  MIN, MAX,  AND  INCREMENTS  IN  WN' 

READ  {*,*)  WNMIN, WNMAX, INCR 
WRITE  (*,*)  'ENTER  WNO' 

READ  ( * , * )  WNO 
WEIGHT=435 . 0 


IZ 

=45.0 

L 

=7.3 

RHO 

=  1.94 

G 

=32.2 

XG 

=0.0104 

MASS 

=WEIGHT/G 

MASS 

=MASS/(0.5*RHO*L**3) 

IZ 

=IZ/(0.5*RHO*L**5) 

XG 

=XG/L 

YRDOT 

=-0.00000 

YVDOT 

=-0.03430 

YR 

=+0.00000 

YV 

=-0.10700 

YDRS 

=+0.01241 

YDRB 

=+0.01241 

NRDOT 

=-0.00047 

NVDOT 

=-0.00000 

NR 

=-0.00390 

NV 

=-0.00000 
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NDRS  =-0 .337*YDRS 
NDRB  =+0 . 283*YDRS 
DH  = ( IZ-NRDOT) * (MASS-YVDOT) - 
&  (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 

All= ( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DH 


A12= ( (IZ-NRDOT) * ( -MASS+YR) - 
&  (MASS*XG-YRDOT) * { -MASS*XG+NR) ) /DH 

A21= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DH 
A22= ( (MASS-YVDOT) * ( -MASS*XG+NR) - 
&  (MASS*XG-NVDOT) * { -MASS+YR) ) /DH 

Bll= ( (IZ-NRDOT) *YDRS- (MASS*XG-YRDOT) *NDRS) /DH 
B12= ( (IZ-NRDOT) *YDRB- (MASS*XG-YRDOT) *NDRB) /DH 
B21= ( (MASS-YVDOT) *NDRS- (MASS*XG-NVDOT) *YDRS) /DH 
B22= ( (MASS-YVDOT) *NDRB- (MASS*XG-NVDOT) *YDRB) /DH 
B1  =B11-B12 
B2  =B21-B22 

Cl= (A11*A22-A21*A12 ) * (A21*B1-A11*B2 ) 

C2= (A11+A22) * (A21*Bl-All*B2 ) +B2* (A11*A22-A21*A12 ) 
C3=- (A21*B1-A11*B2) **2 
A=C1/C2 
B=C3/C2 
C 

EPS=l.D-5 

ILMAX=1500 

C 

DO  1  11*1, INCR 
C 

WN  =WNMIN+  (WNMAX-WNMIN)  *  (II-l)  /  (INCR-1) 

C  print  *,wn 

ALPHAO  =WN*  *  3 
ALPHA1=2 . 15*WN**2 
ALPHA2 =1 . 7  5  *WN 
KPSI=-ALPHA1/B 
KY  = -ALPHAO /B 
KR  =- (ALPHA2+A) /B 
GAMAO  =WNO*  *  3 
GAMA1=2.15*WN0**2 
GAMA2 =1 . 7  5  *WNO 
LY=A+GAMA2 
LPSI=A*LY+GAMA1 
LR=A*LPSI +GAMA0 
C234567 

DO  2  J=l, IC 

CCRI T=CMIN+ (J-l) * (CMAX-CMIN) / (IC-1) 

C  A  MATRIX 

AMAT ( 1 , 1 ) =0 . 0 
AMAT(1, 2 ) =1 . 0 
AMAT ( 1 , 3 ) =0 . 0 
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AMAT(1, 4) =0 . 0 

AMAT ( 1 , 5 ) =0 . 0 

AMAT ( 1 , 6 ) =0 . 0 

AMAT (2, 1)  =0 . 0 

AMAT (2,2) =A 

AMAT(2 , 3 ) =0 . 0 

AMAT (2,4) =B*CCRIT*KPSI 

AMAT (2 , 5) =B*KR 

AMAT (2,6) =B*KY 

AMATO, 1)=1.0 

AMATO, 2)=0.0 

AMATO,  3  )  =0 . 0 

AMATO,  4)=0.0 

AMATO, 5)=0.0 

AMATO,  6)=0.0 

AMAT(4, 1)  =0 .0 

AMAT(4, 2 ) =0 . 0 

AMAT (4 , 3 ) =LPSI 

AMAT (4, 4) =0 . 0 

AMAT (4, 5)  =1 . 0 

AMAT(4, 6) =-LPSI 

AMAT(5, 1)  =0 . 0 

AMAT(5,2)=0.0 

AMAT  (5,3)  =LR 

AMAT (5,4) =B*CCRIT*KPSI 

AMAT (5,5) =B*KR 

AMAT (5,6) =-LR+B*KY 

AMAT  ( 6 , 1 )  =0 . 0 

AMAT (6, 2)  =0 .0 

AMAT (6,3) =LY 

AMAT ( 6 , 4 ) =1 . 0 

AMAT(6, 5)  =0 .0 

AMAT  (6,6)  =-LY 

C 

CALL  RG (6,6, AMAT, WR, WI , 0 , ZZZ, IV1 , FV1 , IERR) 
CALL  DSTABL ( DEOS , WR , WI , FREQ ) 

U=CCRIT 

IF  (J.GT.l)  GO  TO  10 

DEOSOO=DEOS 

UOO=U 

LL=0 

GO  TO  2 

10  DEOSNN=DEOS 
UNN=U 

PR=DEOSNN*DEOSOO 

IF  (PR.GT.0.D0)  GO  TO  3 

LL=LL+1 

IF  (LL.GT.3)  STOP  1000 

IL=0 

UO=UOO 
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UN=UNN 
DEOSO=DEOSOO 
DEOSN = DEOSNN 
6  UL=UO 

UR=UN 

DEOSL=DEOSO 
DEOSR=DEOSN 
U= (UL+UR) /2 . DO 
CCRIT=U 
AMATd, 1) =0 . 0 
AMATd, 2)  =1.0 
AMAT (1 , 3 ) =0 . 0 
AMATd,  4)  =0 . 0 
AMAT (1 , 5) =0 . 0 
AMATd, 6)=0.0 
AMAT (2 , 1 ) =0 . 0 
AMAT (2,2) =A 
AMAT (2, 3)  =0.0 
AMAT (2,4) =B*CCRIT*KPSI 
AMAT (2, 5) =B*KR 
AMAT (2,6) =B*KY 
AMATO,  1)  =1.0 
AMATO, 2)=0.0 
AMATO,  3)  =0.0 
AMAT (3 , 4)  =0 . 0 
AMATO, 5)=0.0 
AMATO, 6)=0.0 
AMAT  (4,1)  =0.0 
AMAT (4, 2 )  =0 . 0 
AMAT (4,3) =LPSI 
AMAT (4, 4)  =0 . 0 
AMAT (4, 5)  =1.0 
AMAT (4, 6) =-LPSI 
AMAT (5, 1) =0 . 0 
AMATd,  2 )  =0 . 0 
AMAT  (5,3)  =LR 
AMAT (5,4) =B*CCRIT*KPSI 
AMAT (5,5) =B*KR 
AMAT (5,6) =-LR+B*KY 
AMATd,  1)  =0 . 0 
AMATd,  2)  =0 . 0 
AMAT  (6,3)  =LY 
AMAT{6, 4)  =1 . 0 
AMAT(6, 5)  =0 . 0 
AMAT (6,6) =-LY 
C 

CALL  RG(6, 6, AMAT,WR,WI, 0, ZZZ, IV1,FV1, IERR) 
CALL  DSTABL ( DEOS , WR , WI , FREQ ) 

U=CCRIT 

DEOSM=DEOS 
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O  H*  to 


UM=U 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

UO=UL 

UN=UM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS (UL-UM) 

IF  (DIF.GT.EPS)  GO  TO  6 

U=UM 

GO  TO  4 

5  IF  (PRR.GT.0.D0)  STOP  3200 
UO=UM 
UN=UR 

DEOSO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS (UM-UR) 

IF  (DIF.GT.EPS)  GO  TO  6 
U=UM 

4  LLL=10+LL 

CCRIT=U 

WRITE  ( LLL , * )  CCRIT, WN 
3  UOO=UNN 

DEOSOO=DEOSNN 
CONTINUE 
CONTINUE 

2001  FORMAT  (215) 

END 

SUBROUTINE  DSTABL ( DEOS , WR , WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR ( 6 )  , WI ( 6 ) 

DEOS=-l . 0D+20 
DO  1  1=1,6 

IF  (WR(I) .LT.DEOS)  GO  TO  1 
DEOS=WR(I) 

IJ=I 

1  CONTINUE 

OMEGA=WI ( I J ) 

OMEGA=DABS (OMEGA) 

RETURN 

END 
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APPENDIX  C 

HOPP  BIFURCATION  PROGRAM  FOR  [X,X]  BASIS 


PROGRAM  KKPSI 

HOPF  BIFURCATIONS 
NOMOTO'S  FIRST  ORDER  MODEL 

CALCULATIONS  FOR  K  AND  CCRITICAL  IF  KPSI  CHANGES  W /  C 


234567 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 


& 

& 

& 

& 

& 

& 

& 

& 


DOUBLE  PRECISION  K1 , K2 , K, LPSI , LY, LR, IZ, L, 

MASS ,  NV,  NR,  NVDOT,  NRDOT,  NDRS ,  NDRB,  KPSI ,  KR,  KY ,  K3 , 
L21 , L22 , L23 , L24 , L3 1 , L32 , L33 , L34 , L51 , L52 , L53 , L54 , 
Mil ,  Ml  2  ,M13,M14,M15,M16,M21,M22,M23,M24,M25,M26, 
M31,M32  ,M33  ,M34,M35,M36,M41,M42,M43,M44,M45,M46, 
M51 ,  M52 ,  M53 ,  M54 ,  M55 ,  M56 ,  M61 ,  M62 ,  M63 ,  M64 ,  M65 ,  M66 , 
Nil ,  N12 ,  N13 ,  N14  ,N15,  N16,  N21,  N22 ,  N23 ,  N24,  N25,  N26, 
N31,N32,N33,N34,N35,N36,N41,N42,N43,N44,N45,N46, 
N51 ,  N52 ,  N53 ,  N54 ,  N55 ,  N56 ,  N61 ,  N62 ,  N63 ,  N64 ,  N65 ,  N66 


DIMENSION  AMAT  (6,6)  ,T(6,6)  ,  TINV(6, 6)  ,  FV1  ( 6)  ,  IV1  (6)  ,YYY[6,6) 
DIMENSION  WR  ( 6 )  ,  WI  ( 6 )  ,  TSAVE  (6,6),  TLUD  (6,6)  ,  IVLUD  ( 6 )  ,  SVLUD  ( 6 ) 


'NEW' ) 


OPEN 

(10, FILE='CVWN1 .RES' , 

OPEN 

(11, FILE='AKPSI .MAT' , 

OPEN 

( 12 , FILE= ' IREAL . MAT ' 

OPEN 

( 13 , FILE= ' RVALS .MAT ' 

OPEN 

( 14 , FILE= ' KKKKK . MAT ' 

WEIGHT=435 . 0 

IZ 

=45.0 

L 

=7.3 

RHO 

=1.94 

G 

=32.2 

:g 

=0.0104 

i-lASS 

=WEIGHT/G 

MASS 

=MASS/ (0.5*RHO*L**3) 

IZ 

=IZ/ (0.5*RHO*L**5) 

XG 

=XG/L 

YRDOT 

=-0.00000 

YVDOT 

=-0.03430 

YR 

=+0.00000 

YV 

=-0.10700 

YDRS 

=+0.01241 

YDRB 

=+0.01241 
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NRDOT  =-0.00047 
NVDOT  =-0.00000 
NR  =-0.00390 

NV  =-0.00000 

NDRS  =-0 . 337  *YDRS 
NDRB  =+0.283*YDkS 
DH  = (IZ-NRDOT) * (MASS-YVDOT) - 
&  (MASS*XG-YRDOT)  *  ( MAS S*XG -NVDOT) 

All= ( ( IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DH 
A12= ( (IZ-NRDOT) * ( -MASS+YR) - 
&  (MASS*XG-YRDOT) * ( -MASS*XG+NR) ) /DH 

A21= ( (MASS-YVDOT) *NV- ( MAS S*XG- NVDOT) *YV) /DH 
A22= ( (MASS-YVDOT) * ( -MASS*XG+NR) - 
&  (MASS*XG-NVDOT) *( -MASS+YR) ) /DH 

Bll= ( (IZ-NRDOT) *YDRS- (MASS*XG-YRDOT) *NDRS) /DH 
B12= ( (IZ-NRDOT) *YDRB- (MASS*XG-YRDOT) *NDRB) /DH 
B21= ( (MASS-YVDOT) *NDRS- ( MAS S*XG -NVDOT) *YDRS) /DH 
B22= ( (MASS-YVDOT) *NDRB- ( MAS S*XG- NVDOT) *YDRB) /DH 
B1  =B11-B12 
B2  =B21-B22 
C 

200  WRITE  (*,1004) 

READ  ( * , * )  IWN 
INCR=IWN 

WRITE  (*,*)  'ENTER  OBSERVER  WN' 

READ  (*,*)  WNO 
205  WRITE ( * , 1007 ) 

READ  (*,*)  A3 
C  DO  is  Dsat 

50  WRITE  (*,1006) 

READ  (*,*)  DO 

C1=(A11*A22-A21*A12) * (A21*Bl-All*B2 ) 

C2= (A11+A22 ) * (A21*B1-A11*B2 ) +B2* (A11*A22-A21*A12 ) 
C3=- (A21*B1-A11*B2 ) **2 
A=C1/C2 
B=C3/C2 
A3=0 . 0 


DO  1  11=1, INCR 


READ  (10,*)  CCRIT, WN 
C  print  *,wn 

ALPHA0  =WN*  *  3 
ALPHA1=2 . 15*WN**2 
ALPHA2=1 . 75*WN 
KPSI=-ALPHA1/B 
KY  =-ALPHA0/B 
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KR  =-  (ALPHA2-. A)  /B 

GAMAO=WNO**3 

GAMA1=2 . 15*WNO**2 

GAMA2=1.75*WNO 

LY=A+GAMA2 

LPSI=A*LY+GAMA1 

LR=A*LPSI+GAMAO 

C234567 

C  A  MATRIX 

AMAT ( 1 , 1 ) =0 . 0 
AMAT (1 , 2 ) =1 . 0 
AMAT ( 1, 3 ) =0 . 0 
AMAT ( 1 , 4 ) =0 . 0 
AMATd, 5) =0 . 0 
AMAT(1, 6) =0 . 0 
AMAT ( 2 , 1 ) =0 . 0 
AMAT (2,2) =A 
AMAT (2 , 3 ) =0 . 0 
AMAT (2,4) =B*CCRIT*KPSI 
AMAT (2,5) =B*KR 
AMAT (2  )=B*KY 
AMAT  (3  =1.0 

AMATO,-;  =0.0 
AMAT(3, 3)=0.0 
AMAT(3, 4) =0.0 
AMAT(3, 5)=0.0 
AMAT(3, 6) =0.0 
AMAT ( 4 , 1 ) =0 . 0 
AMAT  (4, 2  )  =0 . 0 
AMAT (4,3) =LPSI 
AMAT (4, 4) =0 . 0 
AMAT (4, 5) =1 . 0 
AMAT (4,6) =-LPSI 
AMATO,  1)  =0.0 
AMAT (5, 2 )  =0 . 0 
AMAT  (5,3)  =LR 
AMAT (5,4) =B*CCRIT*KPSI 
AMAT (5,5) =B*KR 
AMAT (5,6) =-LR+B*KY 
AMAT  ( 6 , 1 )  =0 . 0 
AMAT ( 6 , 2 ) =0 . 0 
AMAT (6,3) =LY 
AMAT (6, 4) =1 . 0 
AMAT (6, 5)  =0 .0 
AMAT (6,6) =-LY 
DO  11  1=1,6 
DO  12  J=l, 6 

ASAVE ( I , J) =AMAT ( I , J) 
12  CONTINUE 

11  CONTINUE 
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CALL  RG ( 6 , 6 , AMAT , WR , WI , 1 , YYY , IV1 , FV1 , IERR ) 
CALL  DSOMEG ( I EV , WR , WI , OMEGA , CHECK ) 

C  WRITE  (*,*)  IEV 

WRITE  (12, *)  (WR (I REAL) , IREAL=1, 6) 

OMEGA 0  = OMEGA 
DO  5  1=1,6 

T(1, 1)  =YYY  (I,  IEV) 

T(I, 2) =-YYY (I, IEV+1) 

5  CONTINUE 

IF (IEV.EQ.1)  GO  TO  13 
IF  ( IEV .  EQ .  2  )  GO  TO  14 
IF (IEV.EQ . 3 )  GO  TO  15 
IF (IEV.EQ.4)  GO  TO  16 
IF (IEV.EQ. 5)  GO  TO  17 
STOP  3004 

13  DO  21  1=1,6 

T ( I , 3 ) =YYY (1,3) 

T ( I , 4 ) =YYY (1,4) 

T ( I , 5 ) =YYY (1,5) 

T(I, 6) =YYY (1,6) 

21  CONTINUE 
GO  TO  30 

14  DO  22  1=1,6 

T ( I , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,4) 

T ( I , 5 ) =YYY (1,5) 

T ( I , 6 ) =YYY (1,6) 

22  CONTINUE 
GO  TO  30 

15  DO  23  1=1,6 

T ( I , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T ( I , 5 ) =YYY (1,5) 

T(I, 6) =YYY (1,6) 

23  CONTINUE 
GO  TO  30 

16  DO  24  1=1,6 

T ( I , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T(I, 5) =YYY (1,3) 

T ( I , 6 ) =YYY (1,6) 

24  CONTINUE 
GO  TO  30 

17  DO  25  1=1,6 

T ( I , 3 ) =YYY (1,1) 

T(I, 4) =YYY (1,2) 

T ( I , 5 ) =YYY (1,3) 

T(I, 6) =YYY (1,4) 

25  CONTINUE 

30  CONTINUE 
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NORMALIZATION  OF  THE  CRITICAL  EIGENVECTOR 
CALL  NORMAL (T) 

INVERT  TRANSFORMATION  MATRIX 

DO  2  1=1,6 
DO  3  J=l, 6 

TINV (I, J) =0 . 0 
TSAVEd,  J)=T(I,  J) 

CONTINUE 

CONTINUE 

CALL  DLUD (6,6, TSAVE , 6 , TLUD , IVLUD ) 

DO  4  1=1,6 

IF  ( IVLUD ( I ) . EQ . 0 )  STOP  3003 
CONTINUE 

CALL  DILU (6,6, TLUD , IVLUD, SVLUD) 

DO  8  1=1,6 
DO  9  J=l, 6 

TINV(I, J) =TLUD(I, J) 

CONTINUE 

CONTINUE 

CHECK  Inv(T)*A*T 
IMULT=1 

IF  (IMULT.EQ.l)  CALL  MULT (TINV, AS AVE, T,A2) 
IF  (IMULT.EQ.O)  STOP 
P=A2 (3,3) 

Q=A2 (4,4) 

R=A2 (5,5) 

S=A2 (6, 6) 

WRITE (21 , * )  P,Q,R,S 

DEFINITION  OF  Nij 

N11=TINV(1, 1) 

N2 1 =TINV (2,1) 

N31=TINV(3 , 1 ) 

N41=TINV(4, 1) 

N51=TINV (5,1) 

N61=TINV (6,1) 

N12=TINV(1, 2) 

N22=TINV(2 , 2 ) 

N32=TINV(3 , 2 ) 

N42=TINV(4, 2) 

N52=TINV(5,2) 

N62=TINV(6, 2) 

N13=TINV (1,3) 
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N23=TINV(2,3) 
N33=TINV(3,3) 
N43=TINV (4,3) 
N53=TINV (5,3) 
N63=TINV(6, 3) 
N14=TINV (1,4) 
N24=TINV{2,4) 
N34=TINV (3,4) 
N44=TINV (4,4) 
N54=TINV{5,4) 
N64=TINV(6, 4) 
N15=TINV(1,  5) 
N25=TINV<2,5) 
N35=TINV(3 , 5) 
N45=TINV{4, 5) 
N55=TINV(5, 5) 
N65=TINV{6, 5) 
N16=TINV(1, 6) 
N26=TINV (2,6) 
N36=TINV (3,6) 
N46=TINV(4,6) 
N56=TINV(5, 6) 
N66=TINV(6, 6) 

DEFINITION  OF  Mij 

M11=T(1, 1) 

M21=T ( 2 , 1 ) 

M31=T ( 3 , 1 ) 
M41=T(4, 1) 
M51=T{5, 1) 
M61=T(6, 1) 
M12=T(1, 2 ) 

M22=T(2 , 2 ) 
M32=T(3, 2) 
M42=T(4, 2) 
M52=T(5, 2 ) 
M62=T(6,2) 
M13=T(1, 3) 

M23=T (2 , 3 ) 

M33=T(3 , 3 ) 
M43=T(4, 3) 
M53=T(5, 3) 
M63=T(6, 3 ) 
M14=T(1,4) 
M24=T(2, 4) 

M34=T ( 3 , 4 ) 
M44=T(4, 4) 
M54=T(5, 4) 
M64=T(6, 4) 


ooonoooo 


M15=T(1, 5) 

M25=T(2, 5) 

M35=T{3 , 5) 

M45=T (4 , 5 ) 

M55=T{5, 5) 

M65=T(6, 5) 

M16=T(1 , 6) 

M26=T(2, 6) 

M36=T ( 3 , 6 ) 

M46=T(4, 6) 

M56=T(*j,  6) 

M66=Ti6, 6) 

Kl=l . /8 . * ( (ALPHA2**3 ) +ALPHA0 ) / (ALPHA2 ) 

K2=3 . *A3- . 5* ( ALPHA2 **2) / ALPHAO 

K3=l . / ( (B**2) * (D0**2 ) ) * (ALPHA2+A) * ( (ALPHA0/ALPHA2 ) + (A**2 ) ) 
K=K1* (K2+K3 ) 

print  *,  wn,k,ccrit 

BL=B/ (3*D0**2) 

C 

C234567890123456789012345678901234567 8901234567890123456789012 345 
L21=A3*M21**3-BL*(CCRIT**3*KPSI**3*M41**3+KR**3*M51**3+ 

&  KY**3*M61**3+ 

&  3*CCRIT**2*KPSI**2*KR*M41**2*M5I+ 

&  3*CCRIT**2*KPSI**2*KY*M41**2*M61+ 

&  3*CCRIT*KPSI*KR**2*M51**2*M41+3*CCRIT*KPSI*KY**2*M61**2*M41+ 
&  3*KR*KY**2*M61**2*M51+3*KR**2*KY*M51**2*M61+ 

&  6  *  CCRIT*  KPS I * KR* KY*M4 1 *M5 1 *M6 1 ) 

L22=A3*M22**3-BL*(CCRIT**3*KPSI**3*M42**3+KR**3*M52**3+ 

&  KY**3*M62**3+ 

&  3*CCRIT**2*KPSI**2*KR*M42**2*M52+ 

&  3*CCRIT**2*KPSI**2*KY*M42**2*M62+ 

&  3*CCRIT*KPSI*KR**2*M52**2*M42+3*CCRIT*KPSI*KY*"2*M62**2*M42+ 
&  3  *KR*KY*  *2  *M62  *  *2  *M52 +3  *KR* *  2  *KY*M52  *  *2  *M62 + 

&  6*CCRIT*KPSI*KR*KY*M42*M52*M62) 

L23=3*A3*M21**2*M22-BL*(3*CCRIT**3*KPSI**3*M41**2*M42+ 

&  3*KR**3*M51**2*M52+ 

&  3*KY**3*M61**2*M62+ 

&  3*CCRIT**2*KPSI**2*KR* (M41**2*M52+2*M41*M42*M51)  + 

&  3*CCRIT**2*KPSI**2*KY* (M41**2*M62+2*M41*M42*M61) + 

&  3*CCRIT*KPSI*KR**2* (M51**2*M42+2*M51*M52*M41)  + 

&  3*CCRIT*KPSI*KY**2MM61**2*M42+2*M61*M62*M41)  + 

&  3*KR*KY**2*(M61**2*M52+2*M61*M62*M51)  + 

&  3*KR**2*KY* <M51**2*M62+2*M51*M52*M61)  + 

&  6*CCRIT*KPSI**KR*KY* (M41*M51*M62+ (M41*M52+M42*M51 ) *M61 ) ) 
L24=3*A3 *M21*M22**2-BLM3*CCRIT**3*KPSI**3*M41*M42**2+ 

&  3*KR**3*M51*M52**2+ 

&  3*KY**3*M61*M62**2+ 
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fh  «■>  fr*  fr*  ff>  x>  x*  *>  x> 


&  3*CCRIT**2*KPSI**2*KIV  'M42**2*M51+2*M41*M42*M52 ) + 

&  3*CCRIT**2*KPSI**2*KY* (M42**2*M61+2*M41*M42*M62 ) + 

&  3*CCRIT*KPSI*KR**2* (M52**2*M41+2*M51*M52*M42) + 

&  3*CCRIT*KPSI*KY**2* (M62**2*M41+2*M61*M62*M42 ) + 

&  3*KR*KY**2* (M62**2*M51+2*M61*M62*M52)  + 

&  3*KR**2*KY* (M52**2*M61+2*M51*M52*M62)+ 

&  6*CCRIT*KPSI*KR*KY* (M42*M52*M61+(M41*M52+M42*M51) *M62) ) 

L31= ( -1/6) *M11**3 
L32= (-1/6) *M12**3 
L33= ( -1/2 ) *M11**2*M12 
L34=(-1/2)*M11*M12**2 

C2345678901234567 8901234567 8 9 012 34 567 89012 345 67 8 9012 3 456789 012345 
L51=-BL*(CCRIT**3*KPSI**3*M41**3+KR**3*M51**3+KY**3*M61**3+ 
3*CCRIT**2*KPSI**2*KR*M41**2*M51+ 
3*CCRIT**2*KPSI**2*KY*M41**2*M61+ 
3*CCRIT*KPSI*KR**2*M51**2*M41+3*CCRIT*KPSI*KY**2*M61**2*M41+ 
3*KR*KY**2*M61**2*M51+3*KR**2*KY*M51**2*M61+ 
6*CCRIT*KPSI*KR*KY*M41*M51*M61) 

L52=-BL* (CCRIT**3*KPSI**3*M42**3+KR**3*M52**3+KY**3*M62**3+ 
3*CCRIT**2*KPSI**2*KR*M42**2*M52+ 

3  *CCRIT* *  2  * KPS I *  *  2  * KY*M42  *  *  2  *M62 + 
3*CCRIT*KPSI*KR**2*M52**2*M42+3*CCRIT*KPSI*KY**2*M62**2*M42+ 
3*KR*KY**2*M62**2*M52+3*KR**2*KY*M52**2*M62+ 
6*CCRIT*KPSI*KR*KY*M42*M52*M62 ) 
L53=-BL*(3*CCRIT**3*KPSI**3*M41**2*M42+3*KR**3*M51**2*M52+ 
3*KY**3*M61**2*M62+ 

3*CCRIT**2*KPSI**2*KR*(M41**2*M52+2*M41*M42*M51)+ 
3*CCRIT**2*KPSI**2*KY* (M41**2*M62+2*M41*M42*M61) + 
3*CCRIT*KPSI*KR**2*(M51**2*M42+2*M51*M52*M41)+ 
3*CCRIT*KPSI*KY**2* {M61**2*M42+2*M61*M62*M41)  + 

3  *KR*KY*  *2  *  (H61  *  *2  *M52+2  *M61  *M62  *M51 )  + 
3*KR**2*KY*(M51**2*M62+2*M51*M52*M61)  + 
6*CCRIT*KPSI*KR*KY*  (M41*M51*M62+  (M41*M52+M42*M51)  *M61)  ) 
L54=-BL*(3*CCRIT**3*KPSI**3*M41*M42**2+3*KR**3*M51*M52**2+ 

&  3*KY**3*M61*M62**2+ 

&  3*CCRIT**2*KPSI**2*KRMM42**2*M51+2*M41*M42*M52)+ 

&  3*CCRIT**2*KPSI**2*KY* (M42**2*M61+2*M41*M42*M62 )  + 

&  3*CCRIT*KPSI*KR**2* (M52**2*M41+2*M51*M52*M42 )  + 

&  3 *CCRIT*KPSI *KY* * 2 *  (M62 * *2 *M41+2 *M61*M62 *M42 )  + 

&  3 *KR*KY**2 *  (M62 *  *2  *M51 +2  +M61  *M62 *M52 )  + 

&  3*KR**2*KYMM52**2*M61+2*M51*M52*M62)  + 

&  6*CCRIT*KPSI*KR*KY*  (M42*M52*M61+  (M41*M52+M42*M51 )  *M62 )  ) 

R11=N12*L21+N13*L31+N15*L51 
R12=N12*L23+N13*L33+N15*L53 
R13=N12*L24+N13*L34+N15*L54 
R14=N12*L22+N13*L32+N15*L52 
R21=N22*L21+N23*L31+N25*L51 
R22=N22*L23+N23*L33+N25*L53 
R23=N22*L24+N23*L34+N25*L54 
R24=N22*L22+N23*L32+N25*L52 
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c 

WRITE  (13,*)  R11,R13,R22,R24 
K= (3*R11+R13+R22+3*R24) /8 
WRITE  (11,2001)  WN, K, CCRIT 
1  CONTINUE 
STOP 

1004  FORMAT  ('  ENTER  NUMBER  OF  DATA') 

1006  FORMAT  ('  ENTER  DELTAS AT. ' ) 

1007  FORMAT  ('ENTER  A3') 

2001  FORMAT  (3E15.5) 

END 

SUBROUTINE  DSOMEG ( IJK , WR , WI , OMEGA , CHECK ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  WR ( 6 ) , WI ( 4 ) 

CHECK=-1 . OD+25 
DO  1  1=1,6 

IF  (WR ( I ) .LT. CHECK)  GO  TO  1 
CHECK=WR ( I ) 

IJ=I 

1  CONTINUE 
OMEGA=DABS (WI (IJ) ) 

IF  (WI(IJ) .GT.0.D0)  IJK=IJ 
IF  (WI(IJ) .LT.0.D0)  IJK=IJ-1 
RETURN 
END 

SUBROUTINE  NORMAL (T) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  T ( 6 , 6 ) , TNOR (6,6) 

CFAC=T(1, 1) **2+T(l, 2) **2 
IF  (DABS(CFAC) .LE. (l.D-10) )  STOP  4001 
TNOR ( 1 , 1 ) =1 . DO 

TNOR (2, 1) = (T(l, 1) *T(2, 1) +T(2,2) *T(1, 2) ) /CFAC 
TNOR (3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2)) /CFAC 
TNOR ( 4 , 1 )  =  ( T ( 1 , 1 ) *T ( 4 , 1 ) +T ( 4 , 2 ) *T ( 1 , 2 ) )/ CFAC 
TNOR { 5 , 1 )  =  { T ( 1 , 1 ) *T { 5 , 1 ) +T (5,2)*T(1,2) )/ CFAC 
TNOR (6,1)=(T(1,1)*T(6,1)+T(6,2)*T(1,2))/ CFAC 
TNOR ( 1 , 2 ) =0 . DO 

TNOR (2,2)=(T(2,2)*T(1,1) -T(2 , 1) *T( 1, 2 ) ) /CFAC 
TNOR (3 ,2)=(T(3,2)*T(1,1) -T(3 , 1) *T(1 , 2 ) ) /CFAC 
TNOR (4 , 2 ) = (T(4 , 2 ) *T( 1, 1 ) -T(4 , 1) *T( 1 , 2 ) ) /CFAC 
TNOR (5, 2)=(T(5,2)*T(1,1)-T(5,1)*T(1,2)) /CFAC 
TNOR ( 6, 2 ) = (T ( 6, 2 ) *T (1, 1 ) -T ( 6 , 1) *T( 1, 2 ) ) /CFAC 
DO  1  1=1,6 
DO  2  J=l, 2 

T (I , J ) =TNOR ( I , J) 

2  CONTINUE 
1  CONTINUE 

RETURN 
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END 

SUBROUTINE  MULT { TINV , A , T , A2 ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  TINV (6, 6) ,A(6,6),T(6,6),A1(6,6),A2(6,6) 
DO  1  1=1,6 
DO  2  J=l, 6 
A1(I, J)=0.D0 
A2{I,J)=0.D0 

2  CONTINUE 

I  CONTINUE 
DO  3  1=1,6 

DO  4  J=l, 6 
DO  5  K=1 , 6 

A1(I, J)=A(I,K) *T(K, J) +A1 (I, J) 

5  CONTINUE 

4  CONTINUE 

3  CONTINUE 
DO  6  1=1,6 

DO  7  J=l, 6 
DO  8  K=1 , 6 

A2  (I ,  J)  =TINV(I ,  K)  *A1  (K,  J)  +A2  (I ,  J) 

8  CONTINUE 

7  CONTINUE 

6  CONTINUE 
DO  11  1=1,6 

C  WRITE  (*,101)  (A(I,J) ,J=1,6) 

II  CONTINUE 
DO  12  1=1,6 

C  WRITE  (*,101)  (T(I,J),J=1,6) 

12  CONTINUE 

DO  10  1=1,6 

C  WRITE  (*,101)  (A2 (I, J) , J=l, 6) 

10  CONTINUE 

C  WRITE  (*,101)  A2 (1, 1) 

RETURN 

101  FORMAT  (4E15.5) 

END 
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